Photoabsorption spectra in the continuum of molecules and atomic clusters 
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We present linear response theories in the continuum capable of describing photoionization spec- 
tra and dynamic polarizabilities of finite systems with no spatial symmetry. Our formulations are 
based on the time-dependent local density approximation with uniform grid representation in the 
three-dimensional Cartesian coordinate. Effects of the continuum are taken into account either with 
a Green's function method or with a complex absorbing potential in a real-time method. The two 
methods are applied to a negatively charged cluster in the spherical jellium model and to some small 
molecules (silane, acetylene and ethylene). 
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I. INTRODUCTION 

Oscillator strength distribution characterizes the opti- 
cal response of atoms and molecules. Advances in mea- 
surements with synchrotron radiation and high resolution 
electron energy loss spectroscopy have enabled us to ob- 
tain oscillator strength distribution of a whole spectral 
region originated from valence electrons The photon 
energy dependences of molecular photoelectron spectra 
have also been measured. They provide useful informa- 
tion about the dynamics of the photoionization processes. 

Theoretical analysis of photoabsorption spectra above 
the first ionization threshold requires continuum elec- 
tronic wave function in a non-spherical multi-center po- 
tential. Several methods have been developed for this 
purpose (|], including the continuum multiple scatter- 
ing method ||, the Schwinger variational method HJH], 
finite- volume variational method || , the linear algebraic 
method |Q, and the K- matrix method ||. The Stieltjes 
moment method does not directly utilize the scattering 
state but extract continuum spectrum from the spectral 
moments which are calculated with a square-integrable 
basis set Q. The method has been extensively applied 
to large systems [[lO 11 . The complex basis method also 
gives continuum spectra by the calculation with a square- 
integrable basis set Jl2[ . 

The purpose of the present article is to present alter- 
native theoretical methods for the continuum spectra of 
molecules. Our methods rely upon the linear response 
theory in the time-dependent local-density approxima- 
tion (TDLDA) and employ a uniform grid representation 
in a three-dimensional Cartesian coordinate. 

The TDLDA (or alternatively called the time- 
dependent density-functional theory) is an extension of 



the static density-functional theory to describe electronic 
dynamics under a time-dependent external field fl3|| . In 
the practical applications, an adiabatic approximation 
is usually assumed: the same exchange-correlation po- 
tential as that in the static case is used for the time- 
dependent problem. The correlation effect is included 
through the dynamical screening which is represented 
by an induced local potential. In the TDLDA, the con- 
tinuum boundary condition was treated with the radial 
Green's function for spherical systems such as atoms |l4|] 
and clusters in the spherical jellium model [|l5|Jl6| . The 
importance of the dynamical screening effect on the con- 
tinuum response has been stressed. The method has 
been extended for molecules with a single-center expan- 
sion technique 17 1. However, the application was limited 
to small, axially symmetric molecules. 

Since the Hamiltonian in the Kohn-Sham theory is al- 
most diagonal in coordinate representation, a grid rep- 
resentation in the coordinate space provides an econom- 
ical description. A uniform grid representation in the 
three-dimensional Cartesian coordinate which has been 
developed by Chelikowsky et.al [^8) provides a conve- 
nient basis for our purpose. The problem is then how 
to incorporate the scattering boundary condition in the 
uniform grid representation. 

Our first method is based on the real-time method 
which one of the present authors has recently devel- 
oped to calculate dynamic polarizability of finite systems 
19 2(]]. In the method, the time-dependent Kohn-Sham 
equation is solved explicitly in real-time as well as in real- 
space with a uniform grid. The dynamic polarizability of 
a whole spectral region is then obtained at once through 
the time-frequency Fourier transformation. An advan- 
tage of this method is that we do not need to handle 
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matrices of very large dimensionality, since only a single 
Slater determinant is evolved in time. A drawback is that 
the single-particle continuum states cannot be treated ex- 
actly because the electrons are confined in a limited size 
of box (model space). A possible way to avoid the diffi- 
culty would be employing a complex absorbing potential 
at the end of the box, which we shall discuss later. 

Our second method is based on the modified Stern- 
heimer method |2lf| . The method has been widely used 
in the linear response calculations |^2|,^3| . ft recasts the 
response problem into solving the static Schrodinger-type 
equation with a source term. The problem is then to solve 
the equation with an appropriate outgoing boundary con- 
dition. Our recipe here is to solve iteratively the equa- 
tion taking into account the boundary condition employ- 
ing a Green's function of a spherical potential, separating 
the self-consistent potential into long range spherical and 
short range non-spherical parts. 

The paper is organized as follows. We present our 
methods in Sec. |ll|. The real-time method with an 
absorptive boundary potential and the Green's function 
method are explained. In Sec. 
methods will be demonstrated. 



[II 



applications of the 



First the spherical jel- 
lium cluster is considered to confirm reliability of our 
methods. We then show calculations of photoabsorption 
spectra of some small molecules and compare them with 
measurements. We give interpretation of the obtained 
spectra. Finally conclusion are drawn in Sec. IV. 



II. LINEAR RESPONSE IN THE CONTINUUM 
OF THREE-DIMENSIONAL REAL SPACE 

A. Real-time method with an absorbing potential 

In this section we first recapitulate the real-time 
method in TDLDA For systems with relatively 

large number of particles, the real-time method is one 
of the most efficient method to calculate the electronic 
excitations in molecules and atomic clusters. 

The TDLDA equations for a spin-independent N- 
electron system are given in terms of the time-dependent 
Hamiltonian h(t) which is a functional of the density 
n(r, t). 

i^- t 4>i{r,t) = h(t)4>i{r,t), »=l,-..,JV/2 (2.1) 
h(t) = -^-V 2 + V ion 
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We adopt % = 1 throughout the present paper. Vi on is 
an electron-ion potential for which we employ a norm- 



conserving pseudopotential [|25| with a separable approx- 
imation pq] , /i xc is an exchange-correlation potential. 

We represent the electronic wave function on a uni- 
form grid points insi de a certain spatial area fj8[ . All 
potentials in Eq. (2.2) are diagonal in this representation 
except for the nonlocal part of V lQn . The second order dif- 
ferential operator in the kinetic energy is approximated 
employing a nine-point formula. For the time evolution, 
we use an algorithm developed in Ref. f27f] which utilize a 



predictor-corrector method. The integration of Eq. (2.1) 
is approximately carried out with a fourth-order Taylor 
expansion. : 



i(t + At) = e 



_ -ih(t+At/2)At 
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For stable iteration, the time step At should be smaller 
than the inverse of the maximum eigenvalue of the Hamil- 
tonian. The energy and the particle-number conserva- 
tions are then well satisfied. This method was successful 
in nuclear physics to investigate the dynamics of nuclear 
reaction [£7],f28| and has been proved to be fruitful for 
electron systems as well [|l9]]24|] . 

We start with the static LDA problem solving the 
ground-state Kohn-Sham equations and determine the 
initial occupied orbitals <pi(r) and density no(r). The 
conjugate gradient method is used to solve the static 
Kohn-Sham equation. Then, an external field is turned 
on instantaneously at t = 0. In this paper, we shall con- 
sider only the dipole response, adopting V ex t = kor v 8{t) 
(y = x,y, z) which produces an initial state for the time 
evolution as 



<Mr,0+) = e 4fcor -0 4 (r), u = x,y,z. 



(2.5) 



We calculate time evolution with this initial condition 
without any external fields. The dynamical polarizabili- 
ties in the time representation is then obtained as 



e 2 f 

*v{t) = —-r- I d 3 rr^Sn(r,t), v = x,y,z, 
ko J 



(2.6) 



where the transition density is simply given by the dif- 
ference from that of the ground-state 



Sn(r, t) = n(r, t) — n (r). 



(2.7) 



Since the dynamical polarizability is usually defined in 
the ene rgy representation, we take the Fourier transform 
of Eq. 



a u (uj)= / dta u {t)e ZUJt ~ rt/2 , u = x,y,z, (2.8) 
Jo 

where we introduce a smoothing parameter T. In order 
to obtain a(u>) in energy resolution of T, we need to cal- 
culate the time evolution up to T > 27r/r. In the next 
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section, we discuss a gas-phase average of the x, y and z 
direction, as 



*(uj) = 1/3 ^2 a v{w). 



(2.9) 



When the energy is in a region where u> + > 0, a 
part of the electron wave functions 0i(r, t) can escape to 
the continuum. This outgoing part of wave function is 
eventually reflected at the edge of the box and produces a 
spurious discrete structure in the photoabsorption spec- 
tra. Therefore, we now introduce an imaginary absorb- 
ing potential VF(r) to suppress the reflection and mimic 
the continuum. Since the imaginary potential should be 
zero in a region where the ground-state density has a 
finite value, we should adopt a spherical box of large ra- 
dius and the imaginary potential is switched on only in a 
outer shell of width Ar. Although the complex potential 
at the edge of the box slightly violates normalization and 
energy conservation, the box must be large enough to 
preserve the conservations if the ti me e volution is carried 
out using an initial state of Eq. ( |2.5| ) with fco = (the 
ground state). 

We adopt an absorptive potential of a linear depen- 
dence on coordinate, which has been discussed in the 
wave-packet method for molecular collisions |Mp(J: 



W(r) 







for < r < R, 



-iW ^ for R < r < R + Ar. 



(2.10) 



The height Wo and width Ar must be carefully chosen 
to prevent the reflection. There has been an argument 
based on the WKB theory to elucidate a condition of no 
reflection @j30|. 

For simplicity, let us take a one-dimensional system. 
The absorbing potential is given by a form similar to Eq. 
( fUOj ): W(x) = for x < and W(x) = -iW x/Ax for 
x > 0. The WKB solution is 



ip(x) = exp(ikx) + Rexp(—ikx) for x < 0, 

X \-V4 



ip(x) = (1 + R)E l/i [E + iW 
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exp <j l J dxk(x) 



for x > 0, 



(2.11) 
(2.12) 

where k(x) = [2m(E + iWqx/Ax)] 1 / 2 and k = fc(0). 
Conditions we require are \R\ 2 <C 1 (no reflection) and 
|t/>(A:e)| <C 1 (complete absorption). A continuity condi- 
tion for a derivative at x = leads to 

(2.13) 



assuming Wq/E <C 1. The absolute value of wave func- 
tion at x — Ax is approximately given by 



|"0(Aa;)| « exp 



-WoAx(-) 
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(2.14) 



If we demand |i?| 2 < 0.001 and |^(Ax)| 2 < 0.01, we have 

20 = < \W \ < —ArVSmE 3 / 2 , (2.15) 

Arv8m 10 

where Ax has been replaced by Ar. This is the condi- 
tion discussed in Ref. |}0| which we shall test with nu- 
merical calculations. Strictly speaking, the conditions, 
|i?| 2 < 0.001 and \ip{Ax)\ 2 < 0.01, lead to a factor 18.4 
instead of 20 in the left hand side and 0.128 instead of 
1/10 in the right hand side. In actual real-time calcula- 
tions, electrons with different energies are emitted simul- 
taneously. Thus, Ar and Wq must be chosen properly ac- 
cording to the energy spectrum of photoabsorption above 
the ionization threshold. 



B. Green's function method with an outgoing 
boundary condition 

Exact treatment of the continuum is possible with the 
use of a Green's function. For spherical systems, the 
Green's function can easily be constructed by making 
a multipole expansion and discretizing the radial coor- 
dinate. The method was first applied to nuclear giant 
resonance s [^l|j32f| and then applied to photoabsorption 
in atoms [ p3|J14| j. In this section we present a method to 
construct a Green's function in the uniform grid repre- 
sentation for a system without any spatial symmetry. 

The linear response theory is formulated most conve- 
niently using a density-density correlation function |34j ]. 
The general formalism with local-density approximation 
is given in the frequency representation |3^,[l4|] . The tran- 
sition density, which corresponds to the Fourier trans- 
form of Eq. (2.7), can be expressed with the use of the 
independent-particle density-density correlation function 



Sn(r,oj)= / d 3 r' Xo {r,r';uj)V sc f{r',uj), 



(2.16) 



where V^ c f is the self-consistent field which is the sum of 
the external field and the dynamical screening (dielectric) 
field: 



V sc f(r,w) = V cxt (r,uj) + 



,3 / SV[n(r)j 
a r 



5n(r') 



(2.17) 



where T^[rt(r)] is a single-partic le p otential of the time- 
independent version of Eq. ( [2.2| ), defined by h = 
— V 2 /2m + V. Although a part of V is non-local, the 
screening field arises from density-dependent parts of V, 
the direct Coulomb and the exchange-correlation poten- 
tials, which are lo cal. T he calculation neglecting the sec- 
ond term of Eq. Q2.17| ) will be called an "independent- 
particle approximation" (IPA) in the next section. 
The xo has a form 
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Xo(r,r';w) = 2^2 

i n 

+ 4>*(r) 



<^(r)0 ro (r') 

<M r ) ~<?i r 

e-i - e m -u — tr) 

0m(r)^(r') 
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-Mr') , (2-18) 



where r\ is an infinitesimal positive parameter. The fac- 
tor two in the right hand side of Eq. (|2 .18) is the spin 



degeneracy of each single-particle state. Here the summa- 
tion with respect to unoccupied states m contains both 
discrete and continuum states. Instead of taking the ex- 
plicit summation, we may use the single-particle Green's 
function defined by 



(r, r'; E) = <r| (E - h[n Q ] ± i^ 1 |r') 



all 

E 



E - e k ±i-q' 



(2-19) 



The sign determines a boundary condition of the Green's 
function; the (+) for a outgoing wave and the (— ) for an 
incomin g wav e. In fact, the summation with respect to m 
in Eq. ( 2.1§| ) can be extended to all states because the 
summation over the occupied states for the first term 
will be canceled by a contribution from the second term. 
Therefore, using Eq. (2.19), we can rewrite (2.18) as 



X o(r, r'; u) = 2 £ <fc(r) { (^(r, r'; - u) 

i 

+ GW(r,r';e i +a;)}^(r'), (2.20) 

where we assume that the occupied states have real wave 
functions. 

To calculate the dipole photoresponse of a system, one 
should take V ex t = r v , v — x,y,z. Then, the dynamic 
polarizability is given by 

a v {ij) = — e 2 J d 3 rr v 5n(r 1 uu) 1 v = x,y,z. (2.21) 

Since the number of spatial grid points is large, it is 
not convenient and even impossible to construct explic- 
itly the response function xo( r , r '; w) and to perform spa- 
tial integration by summing up over grid points in so lving 
the self-consistent equations, Eqs. (2.16) and (2.17). We 
can avoid the difficulty by converting the integral into the 
equivalent differential equation. We introduce a function 

Mr;E,V scl ) = |A'GW(r,r';£)y sc[ (r')^(r'). (2.22) 

The transition density is then expressed as 



5n(r, u>) 
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+ ^(r;e; + w,F scf )}. (2.23) 



When the energy u> is far below the first ionization thresh- 
old, all the outgoing channels are closed and Eq. ( 2.22j) 



is equivalent to the Schrodinger equation of energy E = 
ti ± u> with a source term, 



{E-h[n ]M) = V sci \<f>i), 



(2.24) 



assuming that the function Mth) outside a box area van- 
ishes. The integral in Eq. ( 2.22 ) is thus converted into a 
linear algebraic equation (2.24). This procedure is known 
as the modified Sternhcimer method |21|. However, since 
this cannot describe a correct asymptotic behavior of the 
wave function, it is not applicable when the energy E is 
close to zero. Furthermore, when E > 0, the method is 
incapable to produce the continuum spectra. 

Thus, the rem aining task is to calculate ipi(r; E, V sc f) 
defined by Eq. ( 2.22 ) under an appropriate outgoing 
boundary condition. For this purpose, we employ an 
integral equation for the Green's function. We start 
with a single-particle problem with a spherical poten- 
tial Vb(r). For instance, the Green's function for free 
particles (Vo( r ) = 0) has an analytic expression, 



G< +) (r,r';£) = - 



2tt 



(2.25) 



where k = \ / 2mE [E > 0). For a negative energy E < 0, 
exp(+ik\r— r'|) in Eq. (2.25) is replaced by exp(— k\v— r'|) 
with k = x/—2mE. The single-particle Green's function 
for a non-spherical (even non-local) potential V(r, r') 
then satisfies the following integral equation 

G (+) {r,r';E) = G { +) {r , r' ; E) + J d 3 r" d 3 r" 

4 +) (r, r" ; E)V(r" ,r")G^ (r'" , r'; E), (2.26) 

where V{r, r') = V(r, r') - V Q {r)5 3 (r - r'). The bound- 
ary condition of G^ determines an asymptotic behavior 
of G^ + \ Substituting this into Eq. ( p. 22 ), we obtain a 
multi-linear equations for ipi( r i E, V sc f), 

xl>i{r;E,V BC () 



d 3 r'd 3 r"G { +) (r, r'; E)V(r', r")^(r" ; E, V scf ) 

= J d 3 r'G ( + ) (r,r' ] E)V sci (r')Mr')- (2.27) 

In solving this equation, we need to evaluate the following 
integral 



*(r;£) 



dVG|, +) (r,r';£)/(r'), 



(2.28) 



where f(r) is either / d s r'V(r, r')^(r'; E, V sc t) or 
V sc f(r)(j)i(r). We note that both of them are zero out- 
side the box. We calculate *£(r;E) again by recasting 
Eq. ( p. 28 ) into equivalent differential equation, 



17 ~ 1 _ 2^ V2 + Vo{r) ' ^ * (r; E) = /(r) - (2 - 29) 
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In the discretization, this is a linear algebraic equation 
for grid points inside the box area. However, the solu- 
tion outside the box area is needed to apply the Lapla- 
cian operator and should be prepared by other method. 
We prepare it by a multipole expansion method. Noting 
that Vb(r) is a central potential, Gq + ^ can be expressed 
in terms of the regular and irregular solutions of radial 
Schrodinger equation in the usual way. 

'max ( + )/ 

*frm r >R = E ' ^ ' Yi m (mim(E), (2.30) 

lm 

$ lm (E) = 2m [ d 3 r' Ul{r '; E) Y lm (r')f(r'), (2.31) 

Jr<R r 

where ui(r; E) and w[ + \r; E) are solutions of the radial 
differential equation being normalized as the Wronskian 
is unity. The ui is regular at the origin and 

w\ + ^ is an outgoing wave at infinity. The l max = 16 is 
chosen in the later calculations. 

We must choose the Vo(r) so as to make the poten- 
tial V"(r) = V(r) — Vb(r) be negligible outside the box. 
For neutral molecules, the attractive ionic potential is ap- 
proximately canceled out by the repulsive direct Coulomb 
term. However we will later employ a gradient-corrected 
exchange potential which possesses a correct asymptotic 
form of — e? jr. Therefore we use a jellium potential, 
ID), with Z = 1 as Vb(r). The w ; (+) (r;^) in Eq. 



( j2.30| ) is an irregular Coulomb wave function with an 
outgoing boundary condition. We use a FORTRAN77 
program "COULCC" in Rcf. [R5| to calculate Coulomb 



functions of complex energies. The ui{r,E) in Eq. (2.31) 



is calculated by integrating the radial Schrodinger equa- 
tion with a fourth-order Runge-Kutta method. 

We now summarize our procedure to solve the response 
equati on. O nce a n exte rnal field I4 x t(r, lo) is given, Equa- 
tions (2.16) and (2.17) constitute a linear equation for 
the transition density Sn(r,uj). Discretizing on a uni- 
form grid, we have a linear algebraic equation with di- 
mensionality equal to the number of grid points. In or- 
der to solve the equation, we need to calculate actions 
of the xo on some functions, which is equivalent to cal- 
culate actions of the G^ + \ ?pi(r; E, Vs C f) is calculated by 
solving another linear equation (2.27). In order to solve 



Eq. (2.27), we need to calculate actions of th e G 



(+) 

This operation is achieved by solving Eq. ( P-29| ), again 
discretizing on grid points. We note that the outgoing- 
boundary condition is imposed at this stage. The wave 
function ^'(r; E) outside the box i s prep ared by a multi- 
pole expansion method, Eqs^J 2^ ) an d (2.31 ). Solutions 
of linear equations ( 2.16[ ), (2.27) and (2.29) are obtained 
by iterative methods. The iterative procedure is summa- 
rized in Fig. ^. 

Our algorithm thus requires to solve multiple nested 
linear algebraic equations. Conjugate gradient (CG) 
method and its variants offer stable and efficient scheme 
to solve these equations. If the energy E is real, Eq. 



( 2.29] ) is a hermitian problem. Therefore, the CG method 
is very powerful to solve the equation. However, if the 
energy is complex, Eq. ( 2.29| ) becomes a non-hermitian 
problem to which the CG method is not applicable. For 
such cases, we adopt a B i-conjugate gradient (Bi-CG) 
method. To solve Eq. (2.27), we have tested different iter- 
ative methods to this non-hermitian problem. We found 
that a generalized conjugate residual (GCR) method is 
one of the most effective algorithms for the present prob- 
lem. For the outer most iteration loop of Eq. ( [2.23| ) , we 
use the GCR method again. 



Self-consistent determination of 8n(r) 
•Iteration for solving 8« = % a (Km+^H 
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Iteration for solving (l-G^tyW = G^Mcr*. 
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■Iteration for solving V 2 /2m- V^¥ -j 
with outgoing boundary conditions 



8n ( " +l, (r) 



FIG. 1. Algorithm to determine the transition den- 
sity Sn(r, ui) is presented. There are three nested iterative 
loops to solve linear equations. 



III. APPLICATIONS 

A. Spherical jellium model for Na»: Test study 

First we have carried out a calculation for a negatively- 
charged Na cluster Naf . We use a spherical uniformly- 
charged jellium potential for Vi ou (r). The jellium poten- 
tial of radius Rj is given by 



V ion (r) 



IRi 



(*)' 



for r < Ri, 
for r > Ri, 



(3.1) 



where Ze is a total charge of a jellium sphere. For the 
Nay , we have Z — 7 and eight valence electrons. The 
main photoabsorption peak is calculated to be above the 
ionization threshold. Thus, this would be a good test 
case to study a response in the continuum and to check 
the applicability of the theories in the previous section. 

We discretize three-dimensional Cartesian coordinates 
with spacing Ax = Ay = Az = 1.5 A and employ grid 



5 



points inside a spherical box of radius R — 12 A. This 
is found enough for the ground-state density to be negli- 
gible at the edge r = R. The number of mesh poin ts is 
2109. The TDLDA Hamiltonian is given by Eq. |Tg ) 
in which the exchange-correlation potential is that of 
Ref. El: 



For the three-dimensional Green's function calculation. 



MxcK r )] 



1.222 
r s (r) 



0.0666 In 1 



11.4 



(3.2) 



n(r). The radius of jellium 



in units of Ry where A/iitr^ 
sphere Ri is 3.93 x 8 1 / 3 a.u. 

Since the jellium potential for Na^~ is spherical, we can 
use the same technique as that of Ref. J3^,[l4| (a Fortran 
program in Ref. Jl6[ | ) and confirm that our methods pro- 
vide the same results. The Green's function is expanded 
in partial waves and is given by 



G| +) (r, r';E) = m (r< ; E)w ( + ] (r> ; E)/W[ui , w^} 



(3.3) 



for the l-th partial wave. Here ui and w\ are regular 
and irregular solutions of the radial differential equation, 
respectively. The boundary condition of u>r at the edge 
r = R is defined as exp(ifcr) with k 2 = 2m(E — e 2 /R) 
for E > e 2 /R and exp(— nr) with n 2 = 2m(e 2 /R — E) 
for E < e 2 / R. Of course, this cause a change of the 
threshold energy by e 2 /R, but it is enough to check va- 
lidity of our methodology. We also add a small imaginary 
parameter Y/2 to the frequency uj, which plays exactly 
the same role as a smoothing para mete r of the Fourier 
transform in real-time method, Eq. (2.8). In this calcu- 
lation, r = 0.1 eV is used and the mesh size for radial 
coordinate is as small as 0.1 A. 




FIG. 2. 



2 3 
co[eV] 

Imaginary (a) and real (b) part of dy- 



namic polarizabilities for Na^T cluster calculated with the 
Green's function method. Solid lines are the results of 
one-dimensional calculation of Ref. [|l(| of jellium sphere 
and square symbols are those of the three-dimensional cal- 
culation without assuming any symmetry. The smoothing 
parameter V = 0.1 eV is used for both calculations. 



in order to justify Eq. (2.26), a condition that V(r) « 
for r > R must be satisfied. We assume that the asymp- 
totic behavior is the same as that of the (radial) one- 
dimensional calculation. Namely, we adopt the Green's 
function of free particles, Eq. ( 2.25| ), with an energy be- 



ing shifted as E — > E — e / R. We also use the same 
smoothing parameter Y = 0.1 eV. It turns out that the 
inclusion of Y helps convergence of iterative procedure. 

We show the results in Fig. ||. The three-dimensional 
calculation has been done for a frequency range < lu < 
5 eV with a step Aui = 0.1 eV. The solid lines are the re- 
sults of one-dimensional calculation in which the Green's 



functions are expanded in the partial wave, Eq. (3.3). 
The squares are the results of the three-dimensional cal- 
culation which perfectly agree with the solid lines. This 
means that the Green's function method on the three- 
dimensional mesh is able to take account properly of the 
continuum effects. 

The static calculation indicates the highest occupied 
molecular orbital (HOMO) —0.37 eV, and a Coulomb 
potential has a value e 2 /R = 1.2 eV at the edge of the 
box R = 12 A. Thus, the ionization threshold is 1.57 eV 
in this calculation. The figure indicates that the pho- 
toabsorption has peaks below and above the threshold. 
In the IPA calculation, we obtain a single peak at an en- 
ergy 1.4 eV. The screening potential brings the peak into 
the continuum region (around 2.35 eV) and splits the sin- 
gle peak into two parts at the threshold. The summed 
oscillator strengths for energies up to 5 eV correspond to 
about 95% of the Thomas-Kuhn-Reiche (TRK) sum rule 
for valence electrons. 




FIG. 3. 



1 2 3 
co[eV] 

Imaginary part of dynamic polarizabilities for 
Naf cluster calculated with the real-time method. The 
absorptive potentials are chosen as (a) Ar = 6 A with 
Wo = 2 eV, (b) Ar = 12 A with Wo = 1 eV, (c) Ar = 21 
A with Wo ~ 1 eV. The smoothing parameter F = 0.1 eV 
is used on the Fourier transform of Eq. (2.5). 
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We now turn to the real-time method. A necessary size 
of space is much larger than that of the Green's function 
method because the absorbing potential is set up in a 
outer region R < r < R + Ar. We prepare t hree kinds 
of absorbing potentials in a linear form of Eq. ( 2.10| ): (a) 
Ar = 6 A with W = 2 eV, (b) Ar = 12 A with W Q = 1 
eV and (c) Ar = 21 A with Wo = 1 eV. Using a square 
mesh of Ax = Ay — Az — 1.5 A, the numbers of mesh 
points are (a) 7153, (b) 17077, and (c) 44473, respec- 
tively. We use a duration of time evolution T = 70 eV -1 
with a time step At = 0.01 eV" 1 and a smoothing param- 
eter F = 0.1 eV for the Fourier transform. Since we have 
a complex potential at the edge of the box, the energy 
is not strictly conserved but the scale of violation turns 
out to be less than 0.1% for this calculation. The imag- 
inary part of dynamic polarizability is shown in Fig. |3j. 
Parts (a), (b) and (c) correspond to the three cases de- 
scribed above. For (a) and (b), the main resonance in 
the continuum region looks like a superposition of three 
different peaks which are located at u> w 2, 2.7 and 3.5 
eV, respectively. However, since we cannot see the same 
behavior in Fig. || (a) , this should be a spurious effect of 
the reflection of outgoing waves. This is confirmed by a 
calculation of (c), for which the structure at oj > 2.5 eV 
is almost disappeared and we see only a smooth tail. 

No w let us check the criterion of a good absorber, Eq. 
( 2.15 ). The position of the main peak is calculated at 
-B res — E'thr = 2.35 — 1.57 w 0.8 eV above the ionization 
threshold. Usi ng E= 0.8 eV and (8m) _1 « 1 eV A 2 for 
electrons, Eq. Q2.15Q reduces to 



^(Ar)" 1 < \W Q \ < 0.07A? 



(3.4) 



where Wq in units of eV and Ar in A. This criterion is 
satisfied for the case (c) but not for (a) and (b). Fig. | 
shows that the real-time method provides us with a cor- 
rect response in the continuum if no reflection occurs at 
the edge of the box. In order to find a suitable absorb- 
ing potential (of a linear form), we find the criterion, 
Eq. ( 2.15| ), very useful even when a condition Wq/E <C 1 
is not satisfied. 



B. Valence shell photoabsorption of silane 

It is well known that the energy of the highest occupied 
orbital does not coincide with the first ionization poten- 
tial in the simplest local-density approximation. This 
fact causes a serious problem for the continuum response 
calculation that the ionization threshold cannot be ade- 
quately described by the static Kohn-Sham Hamiltonian. 
Furthermore, the excited states around the ionization 
threshold appear in too low excitation energies. To rem- 
edy this defect, a gradient correction for the exchange- 
correlation potential has been proposed. We utilize the 
one constructed by van Leeuven and Baerends |$7| which 
we denote as ^ LB \ It is so constructed that the poten- 
tial has a correct — e 2 /r tail asymptotically. The energy 



of the highest occupied orbital also approximately coin- 
cides with the ionization potential. For small molecules, 
TDLDA calculations with this gradient correction have 
shown to give accurate description of discrete excitations 
in small molecules |B8|. 



we em- 



In the following sections III B , QIC and HID, 
ploy a sum of the exchange-correlation potentials of // pz ) 
of Ref. for the local-density part and /i^ LB ^ for the 
gradient correction. We should remark here that an ac- 
curate calculation of the gradient correction jj- 1 ^ (r) be- 
comes difficult at far outside the molecule, because the 
^( LB ) depends on |Vn(r)|/n(r) 4 / 3 which approaches to 
finite but both the numerator and the denominator ap- 
proaches to zero at r — > oo. Thus, we use an explicit 
asymptotic form, ^( LB '(r) = — e 2 /r for r > R c . In the 
followings, the R c is chosen as 6.5 A for silane and 5.5 A 
for acetylene and ethylene. 

When we evaluate the screening field in Eq. ( 2.17 ), 
we neglect the functional derivative of /jA LB )(r) with re- 
spect to the density because this should be a small cor- 
rection. For the real-time calculation, we make the same 
approximation. Namely the time-dependent exchange- 
correlation potential is calculated by 



/ (M .(r./) */ (M .[,;„(r)] - / dV ^ xc(r) 



Mxc[«o(r)] 



5n(r') 
, V PZ) (r) 



5n(r') 



5n(r',t) 



8n(y',t). 



In this section, we discuss the application of our tech- 
niques to a deformed molecule, silane SiEU. We use the 
pseudopotentials which are calculated by the prescrip- 
tions of Refs. |2j|^(|. For Si atom, we employ pseu- 
dopotentials for s,p,d waves and take ci- wave poten- 
tial as a local part. For the hydrogen atom, we take 
s and p potentials with the latter as a local potential. 
The box is a sphere of radius R = 7 A discretized in 
meshes of Aa; = Ay = Az = 0.4 A. Position of the Si 
atom (nucleus) is located at the origin and four hydro- 
gen atoms are at (1.209,0.0,0.855), (-1.209,0.0,0.855), 
0.0, 1.209, -0.855), and (0.0, -1.209, -0.855) in units of 
The symmetry of this molecule belongs to the tetra- 
hedral (T<j) point group. The valence electronic orbitals 
in the ground state are 

(3 ai ) 2 ,(2t 2 ) 6 : l A,. 

In actual calculations, we do not assume a priori the sym- 
metry, but the degeneracies of electron orbitals according 
to Eq. (IIIB) naturally emerge after the minimization of 
the total energy. For this molecule, the dipole response 
does not depend on the direction of an external field. A 
smoothing parameter T = 0.5 eV is used in the follow- 
ings. 

Utilizing the exchange-correlation potential of // pz ) + 



(LB) 



the calculated occupied valence orbitals in the 



ground state are listed in Table |. If we neglect the gra- 
dient correction term /i' LB \ we obtain —13.5 eV for 3ai 
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and —8.3 eV for 2t2, respectively. The photoabsorption 
oscillator strengths for silane calculated by means of the 
Green's function method are shown in Fig. f|. The oscil- 
lator strength df/dui, the photoabsorption cross section 
cr{uS), and the imaginary part of dynamic polarizability 
Ima(w), are related to each other by 



2ir 2 e 2 df 4ttw 



Ima(uj). 



(3.5) 



mc dio 

A sharp peak at 10 eV in the calculation consists of 
bound discrete peaks overlapped by the width T. In the 
experiment jl3 44, ij, we observe a peak at 10.7 eV, the 
width of which is considered to originate from a cou- 
pling of electronic excitation with ionic motion. Since we 
fix the ion coordinates and calculate vertical electronic 
excitations, we cannot describe the width below the ion- 
ization threshold. 

On the other hand, a broad peak at 14.6 eV in the ex- 
periment is above the ionization threshold and the main 
decay channel is an emission of the electron ||,||. The 
calculation well reproduces this peak although the energy 
shifts slightly to the lower (14 eV). Beyond the peak of 
14.6 eV, a group of tiny peaks is observed in the exper- 
iment, which is assigned the Rydberg states with a hole 
(3ai) -1 embedded in the continuum. This fine structure 
is smeared out in the calculation because the smoothing 
parameter T is much larger than the width of each Ryd- 
berg state. The integrated oscillator strengths for ener- 
gies up to 30 eV is calculated to exhaust 87% of the TRK 
sum rule for valence electrons. About 60% of the valence 
shell strengths lie between 10 and 20 eV (Fig. || (b)). 




10 20 
Photon energy to [ eV ] 

FIG. 4. Calculated and experimental photoabsorption 
oscillator strengths of silane. The thick solid line is the 
calculation compared with synchrotron radiation experi- 
ment (thin solid line) |43] and high resolution dipole (e,e) 
experiment (dotted line) . The dashed line is the IPA 
calculation without dynamical screening. The smoothing 
parameter V — 0.5 eV is used in the calculation. Inset: 
An energy region of 13 < id < 15 eV is magnified and 
the total oscillator strengths are decomposed into those 
associated with different occupied valence electrons, 3ai 
(solid line) and 2ti (dashed line). 



The screening field in Eq. ( 2.17 ) significantly reduces 
a dipole polarizability at low frequency. The calculated 
static dipole polarizability is 5.1 A 3 while it is 7.9 A 3 in 
the IPA. The dielectric effects also change the absorption 
spectra. The screening shifts the resonance energies to 
higher energies because the electron-electron correlation 
is repulsive. Roughly speaking, the IPA overestimates 
the absorption spectra at low frequency and underesti- 
mates at high frequency. As you see in Fig. ||, the di- 
electric effects are very important to reproduce the mag- 
nitude and shape of the photoresponse. 

In order to discuss details of each resonance, it is useful 
to calculate a partial oscillator strength of each occupied 
orbital which should be identified with a photoemission 
spectra if the neutral dissociation channel without elec- 
tron emission is negligibly small. As we see in Ref. fl4f| , 
the photoabsorption cross section can be written as a 
sum of contributions from all occupied orbitals. Here we 
define the partial oscillator strengths (df/dui)i as 



du) 
Im 



du> 
2oj 



(3.6) 



., , d 3 rdVVC cf (r)<Mr) x 



{ (<7« (r, r'; (e t - w)*)) * + (r, r'; e, + w)} 
xy s cf(r>i(r'), (3.7) 



where is defined by Eq. ( 2.19 ) except that the sum- 
mation with respect to k is carried out only for unoccu- 
pied states. 




10 20 
Photon energy co [ eV ] 

FIG. 5. (a) Photoabsorption oscillator strengths of 
silane calculated with the real-time method. The thin line 
corresponds to the result obtained with the box of R = 10 
A without a complex absorbing potential, while the thick 
line to the one with a complex potential of Wo = 4 eV 
with Ar = 10 A in addition to a spherical box of R — 7 A. 
In carrying out the Fourier transform, the smoothing pa- 
rameter P = 0.5 eV is used, (b) The integrated oscillator 
strength f(co). The TRK sum rule indicates /(oo) = 8 
for the valence shell photoabsorption. 
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TABLE I. Calculated eigenvalues of occupied valence orbitals and experimental vertical ionization 
potential (IP) in units of eV. The experimental data are taken from Ref. for silane, from Ref. 
for acetylene and from Ref. J42J for ethylene. 





Silane 






Acetylene 






Ethylene 




occ 


cal 


exp 


occ 


cal 


exp 


occ 


cal 


exp 


state 




IP 


state 




IP 


state 




IP 


(3ai) a 


-17.4 


18.2 




-22.4 


23.5 


(2a g y 


-22.8 


23.7 


(2t 2 ) 6 


-12.4 


12.8 


(2a u ) 2 


-18.4 


18.4 


(2&!„) 2 


-18.6 


19.4 








(3<r 9 ) 2 


-16.7 


16.4 


(1M 2 


-16.3 


16.3 








(ln u ) 4 


-12.1 


11.4 


(3a 9 ) 2 


-14.7 


14.9 



(16 39 ) 2 -13.2 13.0 

(1M 2 -H-7 11-0 



Below the ionization threshold, there are three bound 
transitions in our calculations which we find are associ- 
ated with the excitations of 2t2 valence electrons. The 
first two transitions are considered to correspond to the 
shoulders in the measurements p4] , fl5[ and interpreted as 
transitions to the Rydberg states. The broad resonance 
at 14 eV (14.6 eV in experiment) is above the ioniza- 
tion threshold of the 2t2 orbitals. Therefore, the struc- 
ture may be more complex. An inset of Fig. ^ shows 
the partial oscillator strengths of 3ai and 2t2 orbitals 
in a photon-energy region of 13 to 15 eV. The calcula- 
tion suggests that the resonance structure is due to the 
bound-to-bound excitation of 3eti electrons. The exci- 
tation couples with the bound-to-continuum excitations 
of 2t2 electrons through the dynamical screening effect. 
The coupling shifts the excitation energy up by about 1 
eV and brings about the width due to the autoionization 
process. The partial cross section of 2t2 electrons also 
acquires oscillation as a function of the excitation energy 
due to the rapid change of the induced field. 

Now let us examine the applicability of the real-time 
method. The time step is chosen as At = 0.002 eV^ 1 
and the time evolution is calculated up to T = 12 eV -1 . 
We show the results in Fig. || (a). The bound excited 
states is reasonably described in the calculation. On the 
other hand, the calculation without an absorbing poten- 
tial (thin line) apparently provides a wrong response in 
the continuum region. We must adopt an imaginary po- 
tential to remove several spurious resonances. However, 
it is very difficult to treat the ionization energy properly. 
In order to mimic the continuum, the absorber should 
be effective only for a outgoing electron with positive en- 
ergy. The problem is that the static potential has — e 2 /r 
behavior at large r. Electrons with negative energies 
(—e 2 /R < E < 0) may be absorbed by the imaginary 
potential. Therefore, the effective ionization potential 
becomes — ejjoMO — e 2 /R for this calculation. Taking 
R = 7 A, this shift in the ionization potential amounts 
to about 2 eV. The calculated peak in the continu um is 
located at uj — 14 eV. According to the condition ( 2.151) 
with E = 14-12.4 + 2 = 3.6 eV, we adopt the Ar = 10 A 
and Wq = 4 eV. The number of grid points is 321,781 for 



this case. The thick line in Fig. || (a) shows the result. 
The spurious peaks disappear and the result well agrees 
with that of the Green's function method. The oscillator 
strengths near the ionization threshold (11 < u> < 13 eV) 
indicate some discrepancies, which we naturally expect 
from the above argument. In Fig. | (b), integrated os- 
cillator strength is plotted against the energy. Seven out 
of eight (TRK sum rule) unit of strength locates below 
30 eV. 

Finally we should mention the feasibility of computa- 
tion. In order to calculate the absorption spectra of Fig. 
H the real-time method takes about 10 hours using a 
single CPU of a Fujitsu VPP700E. On the other hand, 
the Green's function method takes about 30 minutes for 
each energy. Thus, the real-time method is faster than 
the Green's function method if the response is calculated 
for over 20 different energies. We have carried out the cal- 
culations for 125 frequencies to obtain the smooth line of 
Fig. | 



C. Valence shell photoabsorption of acetylene 

The acetylene molecule, C2H2, has a symmetry config- 
uration of Dooh- This high symmetry has made possible a 
calculation of the Green's function using a single-center 
expansion pj]. Even so, only two kinds of molecular 
orbitals, 3c g and l7r u , which are primarily derived from 
atomic p states have been taken into account in Ref. jlTj , 
because it was difficult to describe the s-derived states in 
the single-center formulation. In the present paper, we 
consider all valence orbitals, including the 2a g and 2a u 
in addition to the above p-derived orbitals, to calculate 
the photoresponses. 

The spherical box is taken as R = 6 A with meshes of 
Ax = Ay = Az = 0.3 A. All atoms are located on the 
z-axis at ±0.601 A for carbon and ±1.663 A for hydro- 
gen. There are ten valence electrons and the calculated 
energies of occupied orbitals are listed in Table |[ Using 
the Green's function method, we obtain the photoabsorp- 
tion oscillator strengths as a function of photon energy, 
shown in Fig. |(| The calculation indicates a sharp bound 
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resonance at uj = 9.6 eV and a broad structure around 
15 eV which seems to be a superposition of three reso- 
nances. The resonance at 9.6 eV strongly responds to 
a dipole field parallel to the molecular axis. The large 
oscillator strengths in the IPA at uj = 5 ~ 8 eV and at 
uj = 12.5 ~ 13.5 eV are shifted to higher energies by 
the dielectric effects. The agreements with the experi- 
mental data are significantly improved by the inclusion 
of the dynamical screening. The static dipole polarizabil- 
ity is also affected significantly: In the IPA calculation, 
the polarizabilities parallel (a||) and perpendicular (ot±) 
to the molecular axis are an = 10.7 A 3 and a±_ = 3.87 
A 3 . The dynamic screening reduces these values to 4.79 
A 3 and 2.77 A 3 , respectively, which well agree with the 
experimental values, aii = 4.73 and a± = 2.87 A 3 @. 

We find some disagreements between the calculation 
and the experiment in Fig. ^. We observe two distinct 
peaks in the experiments for the broad resonance around 
15 eV. However, the calculation indicates three peaks. 
The lowest (13.2 eV) and the highest ones (15.9 eV) are 
related with responses to a dipole field parallel to the 
molecular axis, while the middle one (14.3 eV) is a re- 
sponse to the perpendicular field. We plot the partial os- 
cillator strengths in an inset of Fig. |[ The lowest peak at 
13.2 eV consists of the transition of 3a g valence orbitals. 
The middle peak at 14.3 eV turns out to consist of contri- 
butions of 2a u and 1tt u orbitals. We need an energy shift 
of these strengths by about 1 eV to reproduce the experi- 
ments. An accurate configuration interaction calculation 
with the Schwinger variational method is available for 
this molecule p|| . This calculation succeeds to reproduce 
the double peak structure. The assignments are consis- 
tent to ours: 3a g for lower and 2a u for higher transitions. 



Photoabsorption oscillator strength 

\ total 




20 30 40 

Photon energy to [ eV ] 

FIG. 6. Calculated and experimental photoabsorp- 
tion oscillator strengths of acetylene. See the caption 
of Fig. ^| The experimental data are taken from Ref. 
^l],[l£j|. Inset: An energy region of 12.5 < uj < 17.5 eV 
is magnified and the total oscillator strength (thick line) 
is decomposed into partial oscillator strengths associated 
with different occupied valence orbital, 2a u (solid line), 
3a g (dashed) and 1tt u (dotted). The contributions of 2a g 
electrons are negligible. 
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FIG. 7. (a) Photoabsorption oscillator strengths of 
acetylene calculated with the real-time method. The thin 
line corresponds to the result obtained with the box of 
R = 10 A without a complex absorbing potential, while 
the thick line to the one with a complex potential of 
Wo = 4 eV and Ar = 10 A in addition to a spherical 
box of R = 6 A. In carrying out the Fourier transform, 
the smoothing parameter V = 0.5 eV is used, (b) The 
integrated oscillator strength f(ui). The TRK sum rule 
indicates f(oo) — 10 for the valence shell photoabsorp- 
tion. 

The experiments also indicate discrete Rydberg series 
around 10 eV. We do not see this structure in the cal- 
culation because we have used a smoothing parameter 
r = 0.5 eV which is much larger than the experimental 
energy resolution. In principle, we can calculate these 
Rydberg states since the potential has a — e 2 /r behavior 
at large distance. In order to obtain these fine struc- 
tures, we have to calculate the response with T « and 
the small frequency mesh Auj (order of meV) must be 
adopted. This is beyond a scope of this paper. 

The real-time calculation is carried out with the same 
imaginary potential as we have used for silane (Wo = 4 
eV and Ar = 10 A). We have 635,371 grid points in this 
space. The time evolution is calculated up to T = 12 
eV -1 with a time step of At = 0.001 eV" 1 . The results 
are shown in Fig. ^ (a). Again the absorbing potential 
is an essential ingredient to obtain sensible results in the 
continuum energy region. The calculated photoabsorp- 
tion spectra well agree with those of the Green's function 
calculation except for minor oscillatory behaviors at high 
energy. This spurious oscillations start to appear at an 
energy u> ~ 22 eV . In fact, the condition of a good ab- 
sorber, Eq. ( [2.15 ), is not satisfied in this energy region. 
Namely the potential (Wo = 4 eV) is too week to absorb 
such high energy particles. In Fig. [?] (b), we show the 
integrated oscillator strengths up to 40 eV. One-fourth of 
the TRK sum rule value of valence electrons is in an en- 
ergy region over 40 eV. It is worth noting that the CPU 
time of the real-time calculation is less than one-fifth of 
that of the Green's function calculation. 
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D. Valence shell photoabsorption of ethylene 

The ethylene C2H4 is the simplest organic 7r-system, 
which possesses the D2/1 symmetry. The two carbon 
atoms are on the x-axis at x = ±0.6695 A and four 
hydrogen atoms are in the x — y plane at (x, y) = 
(1.2342,0.9288), (1.2342,-0.9288), (-1.2342,0.9288), 
(-1.2342,-0.9288) in units of A. There are twelve va- 
lence electrons, and calculated eigenenergies of occupied 
valence orbitals are listed in Table Q. 

The calculations are carried out using the same box 
(R = 6 A) and mesh size (0.3 A) as we have used for 
the acetylene molecule. The photoabsorption oscillator 
strengths calculated with the Green's function method 
are shown in Fig. M. The agreement with an experiment 
f49| is excellent. Almost all the main features of pho- 
toabsorption spectra are reproduced in the calculation. 
The observed bound excited states show different pho- 
toresponses according to the direction of the dipolc field. 
The lowest peak at w = 7.6 eV is mainly a response to a 
dipole field parallel to the molecular (C— C) axis. This is 
associated with the excitations of the HOMO 163^ elec- 
trons. On the other hand, states at 9.8 eV respond almost 
equally to a dipole field of x, y, and z direction, to which 
both l&3 ff and 1&3 U occupied orbitals contribute. A small 
peak at w = 11.4 eV is calculated as a resonance of 3a g 
orbital, which may correspond to a small shoulder in the 
experiment. Beyond 11.7 eV, the HOMO electrons are 
in the continuum. The first prominent peak at 12.4 eV 
is a resonance with respect to a dipole field of y direction 
which is in the molecular plane and perpendicular to the 
C— C bond. The excitations of the \bi u and 163^ elec- 
trons are main components of this resonance. In a region 
of 13.2 < uj < 20 eV, the 163^ electrons can be excited 
into the continuum and produce the smooth background 
of oscillator strengths (0.1 ~ 0.2 eV" 1 ). A peak struc- 
ture at 14.6 eV is made of the excitation of l&2u electrons. 
The peak at 16.4 eV is constructed by excitations from 
the 2&i„ occupied orbitals, while the experiment indicates 
the resonance at 17.1 eV. 

The above analysis is consistent with that in the liter- 
ature Q, except for the 12.4 eV peak (11.9 eV in the 
experiment). Our analysis suggests transition of lbs g and 
1&2 U electrons while Ref. [^9| indicates 3a g . 

Comparing the TDLDA calculation with that of the 
IPA (dashed line) , we see again that the dynamic screen- 
ing effect is very important to reproduce the experimental 
data. This feature of the IPA result is consistent with the 
other calculations without the dynamic screening pb]] . 
The calculated dipole polarizability is a = 4.22 A 3 (5.47, 
3.97, 3.23 A 3 for x, y, z direction, respectively), which 
is in a good agreement with the experimental value, 4.22 
A 3 In the IPA calculation, we obtain a = 7.04 A 3 
(10.4, 5.96, 4.76). 

Results of the real-time calculation are shown in Fig. 
^| (a). The box and imaginary potential we use are the 
same as those for acetylene. We can obtain a sensible re- 



sult if we adopt the absorbing potential. However, again, 
a spurious oscillatory behavior is seen at high energy re- 
gion because the absorbing potential is not strong enough 
to erase those high-energy components. 
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FIG. 8. Calculated and experimental photoabsorption 
oscillator strengths of ethylene. The thick solid line is the 
calculation compared with high resolution dipole (e,e) ex- 
periment (dotted line) The dashed line is the IPA 
calculation without dynamical screening. The smoothing 
parameter V = 0.5 eV is used in the calculation. 
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FIG. 9. Photoabsorption oscillator strengths of ethy- 
lene calculated with the real-time method. See the cap- 
tion of Fig. (?] (a), (b) The integrated oscillator strength 
/(w). The TRK sum rule indicates /(oo) = 12 for the 
valence shell photoabsorption. 



IV. CONCLUSIONS 

We have developed methods based on the TDLDA of 
investigating the photoresponses in the continuum for 
systems with no spatial symmetry: (1) The real-time 
method with an absorbing potential, and (2) the Green's 
function method. These methods allow us to treat the 
photoionization and the dynamical screening effects self- 
consistently. For the real-time method, we have tested 
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imaginary potentials of different ki nds a nd found that 
those satisfying a condition, Eq. ( 2.15| ), are good to 
mimic the continuum effect. However, it is very diffi- 
cult to treat the photoabsorption spectra in vicinity of 
the ionization threshold. There are two reasons: One 
is because the condition, Eq. (2.15), requires a large 
model space to treat such a low-energy emission prop- 
erly. The second is that the ionization threshold is not 
correct when the potential has a 1/r behavior at large dis- 
tance. In addition to this, since the condition is energy 
dependent, it is very difficult to construct a good ab- 
sorber for low-energy and high-energy outgoing electrons 
simultaneously. The advantage of the real-time method is 
a computational feasibility. Utilizing the Fourier trans- 
form, we can obtain the spectra over the whole energy 
region at once. 

The Green's function method possesses a capability of 
an exact treatment of the continuum. Using a Green's 
function of Coulomb asymptotic waves, it is also possi- 
ble to investigate the photoresponse near the threshold. 
The main difficulty of the Green's function treatment is 
a heavy computational task. In order to reduce the CPU 
time, we use a complex energy for the Green's function, 
G(uj + iT/2). We have found that the inclusion of the 
imaginary part T facilitates a convergence of iterative 
procedures in the calculation. The V also plays a role of 
lowering an energy resolution of the calculations. Thus, 
we can choose a value of T depending on the energy res- 
olution required in each problem. We would like to em- 
phasize again that the method is capable of calculating 
response functions of many-electron systems, below, near 
and above the ionization threshold in a unified manner. 

The numerical calculations have been performed for 
a spherical Na^ cluster (test study) and non-spherical 
molecules, S1H4, C2H2, and C2H4. Studies of pho- 
toresponse for deformed molecules including both the 
dynamic screening and the continuum effect are very 
few. An exception is a study of nitrogen and acetylene 
molecules (with a high symmetry -Dooh) by Levine and 
Soven Jl7j] using the single-center expansion technique. 
However, only lir u and 3a g occupied orbitals are included 
in their calculation because of limitation of the single- 
center expansion. Since our calculation has been carried 
out on three-dimensional coordinate meshes, we do not 
need any spatial symmetry including all valence orbitals 
in the calculation. We present the photoabsorption os- 
cillator strengths compared with dipole (e, e) and syn- 
chrotron radiation experiments. The agreement is gen- 
erally excellent. The inclusion of the dynamic screening 
turns out to be essential to reproduce the experiments. 
The IPA calculation overestimates the strengths at low 
energy and underestimates at high energy. 

We are strongly encouraged by the success of our meth- 
ods applied to simple molecules in the present paper. The 
applications to more complex systems may become pos- 
sible in near future. 
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